﻿using System;
using System.Collections.Generic;
using System.Text;

namespace Effects.Wave
{
    /// <summary>
    ///  Class provides way to get Fourier tranformation of given signal.
    /// </summary>
    /// <remarks>
    /// <para>
    ///  FFT algorithm uses common programing technique Devide and conquer, so it only works with input array with length equal to some power of 2.
    ///  In addition it uses twiddle factors (complex roots of unity), that are only dependent on input lenght. Therefore FFT and inverse FFT are not just static methods. 
    ///  To use FFT is necessary create new instance of FFT class with given length of input array. In constructor the twiddle factor are pre-computed and used every time Fft method is called.
    /// </para>
    /// <para>
    ///  Implemented FFT algorithm works purely with complex numbers. Other methods (with real inputs or inverse FFT) are only wrappers for complex-input-FFT method.
    ///  We use popular butterfly method TODO better description of method
    /// </para>
    /// </remarks>
    public class FFT
    {
        /// <summary>
        ///  Expected length of input signal.
        /// </summary>
        private int N;
        /// <summary>
        ///  Pre-computed twiddle factors.
        /// </summary>
        private Complex[] twiddleFactors;

        /// <summary>
        ///  Creates new instance and compute twiddle factors
        /// </summary>
        /// <param name="dataLength"> Expected length of input array </param>
        public FFT(int dataLength)
        {
            N = dataLength;
            twiddleFactors = CreateTwiddles(N);
        }

        /// <summary>
        ///  creates twiddle factors - roots of unity
        /// </summary>
        /// <param name="length"> length of input signal </param>
        /// <returns> array of complex numbers containing twiddle factors </returns>
        private static Complex[] CreateTwiddles(int length)
        {
            Complex[] ret = new Complex[length / 2];
            for (int i = 0; i < length / 2; i++)
            {
                ret[i] = new Complex(Math.Cos(-2 * Math.PI * i / length), Math.Sin(-2 * Math.PI * i / length));
            }
            return ret;
        }

        /// <summary>
        ///  Nearest power of 2 greater or equal than parameter N.
        /// </summary>
        /// <param name="N"> integer N </param>
        /// <returns> Nearest power of 2 greater or equal than parameter N </returns>
        /// <remarks>
        ///  <para>
        ///   This method is used to adjust length of input
        ///  </para>
        /// </remarks>
        public static int nexPowerOf2(int N)
        {
            double log = Math.Log(N, 2);
            return (int)Math.Pow(2, Math.Ceiling(log));
        }

        #region FFT
        // methods for FFT

        /// <summary>
        ///  FFT algorithm for real signal array. FFT is computed in place.
        /// </summary>
        /// <param name="indata"> input signal to transform </param>
        /// <remarks>
        ///  FFT algorithm as it is uses purely complex numbers. This method is only wrapper for complex method. 
        ///  It loads the complex parts with zeros, runs FFT and construct result from complex result. 
        ///  Result is 2 times longer than input signal beacause complex parts are interleaved with real parts. In other words array that 
        ///  comes out of the method looks like {real part 0, imag part 0, real part 1, imag part 1, ..... }
        /// </remarks>
        public void Fft(ref double[] indata)
        {
            double[] data;
            // testing indata length
            int N = indata.Length;
            int powerOf2 = nexPowerOf2(N);
            if (N != powerOf2)
            {
                data=new double[powerOf2];
                Array.Copy(indata, data, N);
            }
            else
                data=indata;

            // creating complex data
            Complex[] complexData = new Complex[powerOf2];
            for (int i = 0; i < powerOf2; i++)
            {
                complexData[i] = new Complex(data[i], 0);
            }

            // complex FFT
            Fft(ref complexData);

            //creating real data back from complex result
            indata = new double[2*data.Length];
            for (int i = 0; i < data.Length; i++)
            {
                indata[2 * i] = complexData[i].re;
                indata[2 * i + 1] = complexData[i].im;
            }
        }

        /// <summary>
        ///  FFT algorithm for complex signal array. It uses butterfly method.
        /// </summary>
        /// <param name="data"> Input signal to transform. </param>
        public void Fft(ref Complex[] data)
        {
            int i, j, k, n1, n2,a,m;
            Complex tmp,twdFactor;

            // bit reversal
            j = 0;
            n2 = N/2;
            for (i=1; i < N - 1; i++) 
            {
                n1 = n2;
                while ( j >= n1 ) 
                {
                    j = j - n1;
                    n1 = n1/2;
                }
                j = j + n1;
             
                if (i < j) 
                {
                    tmp = data[i];
                    data[i] = data[j];
                    data[j] = tmp;
               }
            }

            //FFT
            n1 = 0;
            n2 = 1;
            m = Convert.ToInt32(Math.Log(N, 2));
   
            for (i=0; i < m; i++) 
            {
                n1 = n2;
                n2 = n2 + n2;
                a = 0;
     
                for (j=0; j < n1; j++) 
                {
                     twdFactor = twiddleFactors[a];
                     a +=  1 << (m-i-1);
             
                     for (k=j; k < N; k=k+n2) 
                     {
                        tmp = data[k + n1] * twdFactor;
                        data[k+n1] = data[k] - tmp;
                        data[k] = data[k] + tmp;
                     }
                }
            }
        }

        #endregion

        #region inverseFFT
        // methods fro inverse FFT

        /// <summary>
        ///  Inverse FFT algorithm for real signal array.
         /// </summary>
        /// <param name="indata"> Real input signal array to tranform. </param>
        /// <remarks> 
        ///  FFT algorithm as it is uses purely complex numbers. This method is only wrapper for complex method. 
        ///  It assumes that input array contains interleaved real and imaginary parts, In other words it looks like {real part 0, imag part 0, real part 1, imag part 1, ..... }
        ///  It creates complex array and runs complex inverse FFT. As a result it takes only real parts of complex result (therefore result array is 2-times shorter)
        /// </remarks>
        public void iFft(ref double[] indata)
        {
            /* testing data length
             * just throwing exception, to produce something useful, input must be result of FFT
             * and therefore length should be power of 2 */
            int N = indata.Length;
            int powerOf2 = nexPowerOf2(N);
            if (N != powerOf2)
            {
                throw new ArgumentException("Input array - no power of 2 length");
            }
            
            // zero negative frequencies
            for (int i = N / 2 + 2; i < N / 2; i++)
            {
                indata[i] = 0;
            }
            
            // creating complex data
            Complex[] complexData = new Complex[N/2];
            for (int i = 0; i < N/2; i++)
            {
                complexData[i] = new Complex(indata[2*i], indata[2*i+1]);
            }

            // complex iFFT            
            iFft(ref complexData);

            // returning only real parts
            indata = new double[N/2];
            for (int i = 0; i < N/2; i++)
            {
                indata[i] = complexData[i].re;
            }
        }
        
        /// <summary>
        ///  Inverse FFT algorithm for complex singal array.
        /// </summary>
        /// <param name="data"> Complex input signal to transform </param>
        public void iFft(ref Complex[] data)
        {
            // creating conjugate complex numbers
            for (int i = 0; i < data.Length; i++)
                data[i] = data[i].conjugate();

            // FFT
            Fft(ref data);

            // creating conjugates again
            for (int i = 0; i < data.Length; i++)
                data[i] = data[i].conjugate();
            
            // multiply by 1/N
            Complex slashN = new Complex((double)1 / data.Length, 0);
            for (int i = 0; i < data.Length; i++)
                data[i] *= slashN;

        }

        #endregion
    }

    #region HannWindow

    /// <summary>
    ///  Class containing methods fo applying Hann window function.
    /// </summary>
    /// <remarks>
    /// </remarks>
    public class Hann
    {
        /// <summary>
        ///  Applies Hann window function on given signal.
        /// </summary>
        /// <param name="data"> input signal</param>
        /// <returns> result of Hann window applied on data </returns>
        public static double[] HannWindow(double[] data)
        {
            int length=data.Length;
            double[] outData = new double[length];
            for (int i = 1; i < length+1; i++)
            {
                outData[i-1] = (-0.5 * Math.Cos((2*Math.PI*i)/(length+1)) + 0.5) * data[i-1];
            }
            return outData;
        }

        /// <summary>
        ///  Method applies inverse Hann window on given signal.
        /// </summary>
        /// <param name="data"> input signal</param>
        /// <returns> result of inverse Hann window applied on data </returns>
        public static double[] iHannWindow(double[] data)
        {
            int length = data.Length;
            double[] outData=new double[length];
            for(uint i=1;i<length+1;i++)
            {
                outData[i-1] = data[i-1] / (-0.5 * Math.Cos((2 * Math.PI * i) / (length+1)) + 0.5);
                if (double.IsNaN(outData[i-1]))
                    outData[i-1] = 1;
                
            }
            return outData;
        }
    }

    #endregion


}